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for  solitary  wave  propagation,  where  model  results  are  compared  to  exact  results.  It  is 
found  that  the  model  results  are  much  more  accurate  for  high  solitary  waves  than  are 
earlier,  Boussinesq-type  theories,  and  give  good  results  for  waves  so  high  that  they  are 
almost  breaking.  The  capability  of  the  model  to  treat  a  variety  of  situations  is 
demonstrated  for  colliding  solitary  waves,  nonlinear  dispersive  wave  trains,  waves  in 
channels  of  varying  breadth,  and  undular  bores.  Formally,  the  model  incorporates 
nonlinear  long  wave  theory  exactly,  incorporates  enough  dispersion  to  describe  linear 
waves  with  fourth-order  precision,  so  that  both  shallow  water  waves  and  deep  water 
waves  are  included,  and  describes  accurately  waves  for  which  dispersive  and  nonlinear 
effects  are  both  important. 

The  factors  that  have  made  this  development  possible  include  use  of  a  new  theoretical 
conservation-of-velocity  law  and  associated  formalism  that  reduces  the  dimensionality 
of  the  calculations  by  one,  while  retaining  nonlinaer  and  dispersive  effects  to  high  order.  - 
Accurate  simulations  of  high  waves  can  now  be  carried  out  roughly  two  orders  of 
magnitude  longer  for  the  same  computational  investment  than  is  possible  with  previous 
methods  of  comparable  accuracy.  Conversely,  the  model  is  considerably  more  accurate 
and  able  to  treat  a  wider  range  of  physical  conditions  than  models  that  run  as  fast. 


CONTENTS 


1.  INTRODUCTION  .  1 

2.  ANALYSIS .  6 

3.  THE  NUMERICAL  METHOD  . 19 

4.  TESTS  OF  THE  ACCURACY  AND  STABILITY  OF  THE 

CALCULATIONS  .  26 

5.  EXAMPLES  .  32 

6.  DISCUSSION  .  59 

REFERENCES  .  61 


J 


A  UNIFIED  MODEL  FOR  THE  EVOLUTION  OF  NONLINEAR  WATER  WAVES 


1.  Introduction 

This  paper  describes  a  method  for  obtaining  accurate  representations  of 
evolving  water  waves  under  a  wide  variety  of  physical  situations.  The  model 
gives  accurate  results  for  waves  almost  to  the  breaking  point,  for  waves  in 
channels  of  varying  section,  and  for  waves  throughout  the  wavelength  regime 
from  shallow  water  waves  to  essentially  deep  water  waves.  The  developments 
that  have  made  this  possible  include  the  use  of  a  new  theoretical 
conservation  law  and  associated  formalism  that  reduces  the  dimensionality  of 
the  caculations  by  one,  while  still  retaining  nonlinear  and  dispersive 
effects  to  high  order.  Accurate  simulations  of  high  waves  can  now  be 
carried  out  roughly  two  orders  of  magnitude  longer  for  the  same 
computational  investment  than  is  possible  with  previous  methods  of 
comparable  accuracy  (Fritts  and  Boris,  1979;  barlow  and  Amsden,  1971;  Chan 
and  Street,  1970).  Conversely,  the  model  is  qualitatively  more  accurate 
than  models  that  run  as  fast  (Chwang  and  Wu,  1976;  Peregrine,  1966,  1967; 
Madsen  and  Mei,  1969).  Two  distinct  elements  go  into  the  research:  a 
general  mathematical  analysis  of  the  physical  situation,  and  a  partly 
empirical  formulation  of  an  appropriate  numerical  method  to  generate 
solutions  to  initital /boundary  value  problems  of  interest. 
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The  word  "order"  will  occur  frequently  in  the  paper.  Here,  and 
everywhere  else  in  the  paper,  except  where  explicitly  stated  otherwise,  it 
refers  to  an  ordering  of  parameters  that  measure  the  smallness  of  wave 
amplitude  and  the  smallness  of  wave  train  dispersion.  The  lowest  order 
corresponds  to  long  waves  of  infinitesimal  amplitude  that  travel  through 
water  of  constant  depth.  Specifically,  at  lowest  order  the  wave  has  an 
amplitude  that  is  negligible  when  compared  to  the  undisturbed  depth,  and  any 
relevant  horizontal  scale  of  the  wave  (wavelength  or  half-width)  is 
infinitely  greater  than  the  undisturbed  depth.  The  waves  are  non- 
dispersive.  All  disturbances  can  be  described  as  a  superposition  of  waves 
of  unchanging  shape.  The  wave  shapes  are  arbitrary. 

Three  bodies  of  theory  have  been  developed  that  can  describe  evolving 
waves  that  are  not  so  small  or  not  so  long.  Linear  wave  theory  treats 
infinitesimal  amplitude  waves  of  arbitrary  length,  and  a  satisfactory 
account  of  the  basic  theory  is  available  in  many  places  (see,  e.g.  Lamb 
1932).  Except  at  the  long  wave  limit,  waves  are  dispersive,  and  so  only- 
sinusoidal  waves  of  infinite  extent  can  be  waves  of  unchanging  shape.  Long 
wave  theory  treats  waves  infinitely  long  compared  to  the  depth;  they  can 
have  arbitrary  amplitudes.  'Io  wave  can  propagate  with  unchanging  shape 
(rates  of  change  of  shape  tend  to  vanish  in  the  limit  of  small  amplitudes, 
or,  more  precisely,  small  slope).  Waves  of  depression  car.  be  accommodated 
by  the  the  theory  for  indefinitely  long  times,  but  waves  of  elevation 
steepen,  until  at  some  finite  time  the  horizontal  scale  at  a  wave  front  is 
not  particularly  long  coirpared  to  the  depth.  Extra  features  such  as  bores 
(the  water-wave  analog  to  shock  waves  in  gas  dynamics)  are  added  to  the  long 
wave  theory  to  complete  the  description.  Again,  accounts  of  the  basic 
theory  are  widely  available  (see,  e.g..  Stoker  195?). 


A  third  situation  can  exist,  in  which  effects  of  nonlinearity  and 
dispersion  are  balanced  sufficently  that  waves  can  propagate  for  long 
distances  without  significant  change  of  form,  even  in  the  absence  of 
dissipation.  The  first  person  to  report  them  in  scientific  annals  was 
Russell  (1838,131*5).  He  first  saw  a  solitary  wave  in  a  canal,  and  being  on 
a  horse  at  the  time,  was  able  to  chase  it  along  its  tow  path  for  a  mile  or 
more.  He  later  performed  laboratory  experiments  and  reported  an  accurate 
relationship  between  the  wave  speed  and  wave  amplitude.  Much  theory  has 
gone  into  describing  such  waves,  mostly  published  in  two  spurts.  During  the 
latter  part  of  the  19th  century  the  governing  equations  and  some  solutions 
for  fairly  long  waves  that  can  travel  both  ways  (3oussinesq,  1871, 18^2)  and 
one  way  (Korteweg  and  deVries,  1895)  were  set  down.  In  their  pioneering 
numerical  study  of  solutions  to  the  Korteweg  and  deVries  (KdeV)  Equation, 
Zabuskv  and  Kruskal  (1965)  discovered  that  overtaking  solitary  waves  emerged 
with  properties  no  different  from  their  pre-collision  ones  (apart  from 
roundoff  errors  and  slight  phase  shifts).  They  suggested  that  solitary  wave 
solutions  to  the  KdeV  Equation  be  called  "solitons"  to  emphasize  their 
ability  to  survive  nonlinear  collisions.  Much  research  on  solitons  has 
resulted  (see  the  review  by  Zabusky,  19fil,  for  a  modern  summary). 

The  Boussinesq,  KdeV,  and  P.LW  (Regularized  Long  Wave — an  alternative  to 
KdeV)  descriptions  each  include  the  lowest  order  effects  of  dispersion  and 
nonlinearity.  Each  has  waves  of  unchanging  shape — solitary  waves  and 
periodic  waves.  Each  can  describe  evolving  waves.  But  the  accuracies  of 
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each  description  are  limited  by  their  retaining  only  the  lowest  order 


nonlinear  and  dispersive  terras  of  the  governing  equations. 

Some  researchers  have  worked  hard  on  the  problem  of  actual  solitary 
waves  in  water  (restrictions  on  wave  amplitude  are  removed,  dissipation  is 
not  treated,  and  the  waves  are  assumed  to  be  of  unchanging  shape  from  the 
start).  These  theories  produce  results  on  the  structure  of  solitary  and 
periodic  waves  all  the  way  to  breaking,  and  are  excellent  benchmarks  for 
assessing  the  accuracies  of  theories  of  waves  that  evolve  to  solitary  waves 
But  these  are  special  calculations,  limited  to  waves  of  unchanging  shape. 

The  attempt  here  is  to  develop  a  satisfactory  methodology  that  is 
deraonstratably  accurate  for  waves  where  dispersion  and  nonlinearity  are  in 
approximate  balance,  while  still  retaining  the  ability  to  treat  waves  where 
they  are  not  in  balance.  The  analysis  retains  exactness  for  as  long  as  is 
practical.  At  a  crucial  step  the  choice  of  solving  for  a  flow  field  by  a 
low  order  expansion  has  to  be  made,  and  this  limits  the  description  of 
solitary  waves  to  one  order  better  than  the  Boussinesq  and  Korteveg-deVries 
descriptions  (the  accuracy  turns  out  to  exceed  that  expected  at  the  same 
order  by  comparison  with  exact  solitary  wave  results).  Long  wave  theory  is 
retained  (though  some  tricks  are  necessary  to  maintain  stability),  and  the 
linear  theory  is  retained  through  a  high  order  in  a  dispersion  par?.*"®ter, 
and  can  be  extended  to  even  higher  orders  easily,  if  desired.  The  method  i 
straightforward  and  computationally  efficient. 

Section  2  gives  the  analytical  development.  Section  2A  discusses  the 
pair  of  prognostic  equations  used  in  the  work.  One  of  them  comes  from 
recent  work  (Witting  and  McDonald,  1?P2).  These  exact  equations  io  not 
close  the  system,  however,  and  high  order  expansions  are  made  to  close  it 
(Section  2B).  In  Section  2B  the  fluid  velocity  at  the  bottom  boundary  is 


chosen  as  one  of  the  basic  dependent  variables  to  make  the  expansion  far- 
reaching  (for  linear  waves  the  expansions  converge  for  all  wave  numbers). 

In  Section  29  further  transformations  are  made  that  retain  high  accuracy 
while  keeping  small  numbers  of  terms  in  the  expansions  fthis  is  analagous  to 
representing  a  function  by  Fade  approximants) .  Section  3  outlines  the 
numerical  method  used  to  generate  solutions,  and  discusses  boundary  and 
initial  conditions.  The  numerical  method  is  leapfrog,  which  employs 
centered  time  and  space  finite  differences.  Its  stability  properties  are 
discussed.  Section  4  describes  tests  of  the  accuracy  of  the  model  by 
comparing  properties  of  solitary  waves  that  evolve  from  the  calculation  with 
theories  that  describe  nonevolving  solitary  waves  precisely.  The  model 
solitary  waves  turn  out  to  be  suprisingly  accurate  almost  to  breaking 
amplitudes.  Section  5  describes  the  capability  of  the  model  to  treat  a  wi  ie 
variety  of  physical  siutations.  In  (5A)  the  problem  of  a  head-on  collision 
of  solitary  waves  is  examined.  The  solitary  waves  are  found  to  he  near- 
solitons  (but  not  exact-solitons).  In  (53)  dispersive  wave  trains  are 
generated  to  demonstrate  how  the  model  handles  waves  of  varying  wavelengths. 
Section  5C  shows  simulations  of  laboratory  experiment ~  on  the  propagation  of 
waves  through  channels  of  varying  breadth.  The  simulations  mirror  certain 
aspects  of  the  experiments  that  other  theories  miss.  Section  53  shews 
examples  of  undular  bores  running  through  times  and  distances  long  enough  to 
match  experiments  performed  by  Favre  (1935).  Finally,  Section  6  discusses 
the  key  features  of  the  model  that  are  felt  to  be  chiefly  responsible  for 
its  accuracy,  efficiency,  and  capability. 
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2.  Analysis 


A.  The  exact  governing  equations 

The  physical  model  used  here  takes  the  fluid  to  be  incompressible  and 
inviscid,  undergoing  only  irrctational  motions.  Rigid  impermeable  vails 
bound  the  sides  and  bottom  of  the  fluid.  A  constant-pressure  surface  bounds 
its  top.  Ihis  physical  model  is  that  most  frequently  taken  by  other 
researchers.  Figure  1  shows  the  geometry  and  the  definition  of  some  of  the 
physical  variables. 

Three  exact  relations  hold  for  fluid  motions  in  a  channel  having 
breadth  b  and  instantaneous  depth  h: 
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where  the  horizontal  displacement  is  x,  the  time  t,  the  elevation  above 
still  water  level  n,  the  surface  slope  n' ,  the  components  of  surface 
velocity  ug  (horizontal)  and  v  (vertical),  and  vertically  averaged 
horizontal  velocity  u,  and  the  acceleration  of  gravity  g.  The  vertical 
coordinate  y,  positive  upwards,  does  not  appear  explicitly  in  (2.1-2. 3). 

Equation  (2.1)  is  the  equation  of  continuity  (Stoker,  195"  gives  a 
complete  derivation).  Strictly  speaking,  the  product  nb  is  the  area  above 
still  water  level  and  the  product  uhb  is  the  volume  flux  of  fluid.  Thus, 
should  be  interpreted  as  the  average  elevation  above  still  water  level,  an 
u  as  the  average  horizontal  velocity  over  a  cross  section.  rcr  sidewalls 
that  are  close  enough  together  (b  <<  h),  and  possibly  for  wide  but  gently 
varying  b,  n  and  u  are  sufficiently  constant  across  the  section  that  n  can 
be  taken  as  the  surface  elevation  and  u  the  vertically  averaged  velocity 
(see  Breen  1333).  This  paper  so  interprets  n  and  u,  and  ignores  cross¬ 
channel  variations  throughout.  Equation  (2.2)  can  be  derived  from 
Bernoulli's  Law  for  irrotationai  notions  or  directly  from  the  Euler 
Equations  for  inviscid  but  not  necessarily  irrotationai  flows  (Witting  and 


A)  PLAN  VIEW 


B)  SIDE  VIEW 
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C)  DEFINITION  OF  VARIABLES 

v 


figure  1  -  The  geometry  of  the  cornnut atonal  channel.  Variable  depths  as 
shown  in  (b)  are  not  treated  in  this  paper. 
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the  term  3n/3t  in  (2.5)  can  fce  replaced  with  a  3/3x  term,  using  (2.k)  and  s< 
each  is  of  the  form 
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where  f  is  either  n’o  or  q  ,  and  a' x,t'  contains  n,  3  and  u.  Hosing  the 
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system  thus  demands  that  a  relation  be  found  connect  inn  n,  q  and  u.  This 
relation  is  found  by  considering  the  velocity  potential,  which  satisfies 

72?  =  0  (2.” 


where  for  irrotational  flows 


u'x,y,t  =  3»/3x;  v(x,y,t)  =  3$/3v 
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B.  A  series  solution 


In  terns  of  the  flui  i  velocity  at  the  botton  (taken  here  to  lie  along 
y=C)  u^  i  $o'(x,0,t',  (2."')  gives  the  Taylor  series  expansion 
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Consequently  both  u  and  qg  can  be  written  as  function  of  u^.  Equations 

(2.10)  and  (2.1l),  along  with  (2.U)  and  f2. 5),  close  the  system  and  are  an 
exact  representation  for  irrotatior.al  two  dimensional  waves  in  an 
incompressible  inviscid  fluid  lying  atop  a  rigid  impermeable  horizontal  bed 
If  the  channel  walls  and/or  still  water  depth  vary  only  gradually'  over  a 
characteristic  distance  associated  with  the  wave  motion  (wavelength  or 
solitary  wave  thickness,  for  example)  we  expect  that  the  set  will  produce 
results  that  are  exact  in  the  limit  of  gradual  variation  and  may  be  quite 
accurate  for  sufficiently  gradual  variation  ( f-reen  lq3?). 

C.  Approximations 

In  the  numerical  work  (2.M  and  (2.5''  are  used  to  advance  n  and  q  by 

s 

one  time  level  in  a  finite  difference  approximation.  The  problem  then 
arises  of  how  to  evaluate  the  fields  in  the  3/3x  terms  of  (2.1'  and  f'  a 
this  new  time  level.  Using  the  formulation  here,  this  involves  solving 

(2.11)  for  u-^  from  the  new  values  of  qs  and  h  =  hQ  +  n.  The  method  of 
solution  adopted  here  employs  a  fast  "Tridiagonal  'latrix  Solver"  (Boris, 
197*$).  It  solves  (inverts)  matrices  containing  only  on-diagonal  and  the 
adjoining  off-diagonal  elements.  Thus,  the  finite  difference  form  of  (2.11 
can  be  inverted  if  truncated  after  the  term  -1/2  (h2u^')’.  Moreover, 
repeated  operations  with  the  tridiagonal  solver  can  solve  equations  that 
resemble  (2.1l)  and  perfectly  account  for  the  linear  terms. 
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Rather  than  try  to  retrace  the  reasoning  that  led  to  the  formul 
eventually  adopted,  I  simply  state  the  procedure  now  ar.  1  later 
(2.10)  and  (2.1l). 

The  depth  h  that  appears  in  (2.10)  and  (2.11)  is  expressed 
h=h  +n.  Equation  (2.11)  is  then  rewritten: 


q,  =  u  ^  (h.  2u  '  )  '  -  ^  !h2-h2  )u  '  j  •  -t  ... 
3  b  2  o  n  2  oh 


A  new  velocity  variable  u  is  then  introduced  that  ra.nsforrts  (2 
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q  =  u  -  a  '  (h  2  u '  ) '  -  or  I  (h2-h.  2^u  '  ' 

s  m,n  o  2  ’  o' 


ta  (2)  0 


P  ■  •  fv,  4  ,,  »  »  »  N  »  ».  ^ _  (  fv.4  v«  *+  V,  »  i  ?  1  » 

»•  \  n  u  /  +  ^  i  i  \ .  i  /U  —  •  •  • 

o  2^  c 


■1)  (2) 

where  a'  and  a  are  constants  that  enter  the  linear  terms  o 
m,n  ra,r. 

Equation  (2.10)  is  similarly  decomposed  i.e. : 


u  =  u  -  V1'  h2  u"  -7  (h2-h2)  u’  ’  + 
n,n  o  oo 


+  b'2^  h-  uiY  +Ti-  (h^-hM  uiv 
m,n  o  120  o 


at  ion 

connect  it  to 

as 

(2.12) 

.12)  to  th° 


•2. 1 • ' 

f-  '  t  ’  "3  ' 

.  •  —  •  -4.  _  • 
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The  coefficients  a, 
accuracy  of  linear 

1st  order: 


2nd  order: 


3rd  order: 


Uth  order: 


and  b  are  selected  to  maximize  the 
n,n  m,n 

waves.  The  first  few  combinations  are: 


aa)  =  1/3;  J'2}  =  0; 

o,l  o,l 


»;H  •• 


f  i }  ,  (?) 

H,i  =  2/5i  ai,i  =  C; 


bj1]  =  1/15;  a,(2J  =  0 ; 

L,I  L  ,1 


/ 1  \ 

q  '  =  7  /7  ■ 

al,2  "  ’ 


,(1) 

2,2 


=  V9; 


(?) 

a.  '  =  1  '1  rc  . 


( \) 

51,2  =  2/21;  V?  = 


,(2)  =  i  /k3. 
*2,2 


a  (  3  ^  =  o. 


■3)  n 

a.  .  =  0 


*3  T>  ”  ’  *  *  * 

-  1  -* 


(2.15) 


'2.1c) 


/  o  1  t  N 


'2.18' 


'2,2 


f  ; 3) 

=  1/9;  =  -/^5;  Vo 
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The  dispersion  relation  for  linear  waves  becomes  of  the  form: 


(1st  order) 


1  +  ij  k2h2 
3  o 


(2rd  order) 


1  +  k2h2 
IS  o 

1  +  —  k2h‘~ 

5  o 


(3rd  order) 


*  It 


+  -|  k2h2  +  k4h4 
7  o  105 


( 4th  order) 


+  ±  k2h2 
9  o 


+  ^  kh 


D 


1  +  k2h2  +  ftp  k4'-4 

9  o  63  o 


These  are  the  entries  to  a  Fkde  table  P(N,!<)  or  P'T-1,.")  representing 
tanh  (khQ),  and  are  correct  to  the  order  of  k2hQ2  indicated  in  the 
parenthesis. 

To  second  order  (in  dispersion)  the  calculations  can  be  explained  nos* 
simply.  Use  one  tridiagonal  solve  on  (2.13): 


i  ■  u  -  T  (hV  ) '  -  4  1  (h2-h2'u  ’  1  ' , 
s  5  o  2  o 


(2.20'' 


Shus  obtaining  u  =  u  (x.t  ).  Then  substitute  tnis 
=  new 


value  of  u  into  (2.14): 


u  (x.t  )  =  u  -  h2  u  "  -4  ( h 2 -h 2 '  u 

new  15  o  o 


,,  »  * 


;o  evaluate  the  new  value  of  u. 
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r 


If  the  fields  r  and  q  ,  are  regarded  as  arising  from  the  dual 
expansion  in  an  amplitude  parameter  £=n  /h 0  and  a  dispersion 
parameter  u=h 2/i2,  where  l  is  a  character  1st ic  horizontal  scale  of  th.e  wave 
motion,  then  the  procedure  using  (2.20)  and  retains  all  terms  of  the 

exact  equations  (?.ul  and  (2.5)  that  are  the  following  orders: 

(governs  waves  in  long  wave  theory; 
if  only  terms  in  ey1 22  are  retained, 
we  have  linear  long  wave  theory) 

(an  approximation  to  linear  waves' 

e2y1/2and  ey3/2  (the  Boussinesq  Equations  retain 

terms  only  through  these  orders) 

e3u1/2;  e2y3/2;  ey572  (the  first  order  beyond  the 

Boussinesq  Equations' 


.1/2 


eu 


5/2 


The  lowest  order  terms  that  are  dropped  are: 


0(ey7/2);  0(  e2u3  22  ' ;  2(e3u3',2'  (2.22' 

”oving  to  a  3rd  and  1th  order  for  mi  1 at i o  n ,  which  involves  solving 
(2.13)  with  2  successing  tridiagonal  matrix  solvers,  pushes  the  dispersive 
errors  from  0(sy7/2)  to  0(eu3/2)  and  0(ey11/2)  respectively'. 
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Unfortunately,  the  other  terns  of  '2.22) 


i.e 


,  0  (  £2  an  i  £^u^/'“) 


renain.  Even  so,  the  formulation  stands  at  one  order  in  £  ar.i  u  further 
along  than  the  3oussinesq  Equations.  As  we  shall  see,  the  numerical  result 
yield  solitary  waves  of  a  practical  accuracy  beyond  that  expected  from 
(2.22). 

Figure  2  displays  the  connection  between  this  theory  and  other  water 
wave  theories.  As  the  figure  shows,  each  theory  is  valid  over  a  restricted 
range  of  nonlinearity ,  dispersion,  or  both.  The  merit  of  this  theory  is 
two-fold:  first,  it  incorporates  all  of  long  wave  theory,  no  natter  how 
nonlinear,  and  goes  all  the  way  to  the  fourth  order  in  fc2h  2  in  a  long 
wave  expansion  of  linear  theory.  Second,  this  theory  goes  one  order  beyond 
the  theories  of  3oussinesq,  Korteweg  and  deVries,  and  other  investigators 
who  consider  nonlinear  dispersive  waves.  Consequently,  the  formal  (and,  as 
will  be  shown,  the  practical)  accuracy  of  this  theory  is  higher  than  the 
3oussinesq  type  theories,  and  can  be  applied  to  higher  waves.  In  sunnary, 
this  theory  can  be  called  a  "unified"  theory  of  water  wave  propagation  in 
that  it  incorporates  long  wave  theory,  a  high-order  version  o-p  linear  wave 
theory,  and  can  treat  higher  waves  than  any  other  time^ependent  theory  tha 
includes  both  nonlinearities  and  dispersion. 

To  show  how  well  the  theory  deals  with  linear  waves  that  are  no4-  so 
long.  Figure  3  gives  dispersion  relations  from  various  theories.  The  "2nd 
order"  theory  is  shown  to  be  fairly  accurate.  The  "uth  order"  theory  that 
makes  use  of  (2.1 3)  can  hardly  be  distinguished  from  exact  linear  theory" 
over  the  range  of  kh_  displayed.  This  range  to  kh  =8  takes  us  into 
essentially  deep  water  waves.  For  example,  the  phase  speed  of  waves  having 
khQ=3  is  within  C.0CC0002  of  the  speed  of  waves  in  infinitely"  deep  water. 
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LEGEND 

=3  Li  NEAR  WAVE  THEORY 
CED  LONG  WAVE  THEORY 
SS  BOL'SSINESQ  THEORIES 
GZ3  THIS  THEORY 


0  12  3  4  5  6 

DISPERSIVE  index,  m 


Figure  2  -  Various  water  wave  theories  for  a  rectangular  channel.  Theories 
of  evolving  waves  in  water  generally  are  limited  by  truncating  terms  in  the 
exact  governing  equations  that  involve  dispersion  or  nonlinearity  or  both. 
The  convention  used  here  is  to  consider  the  linear  long  wave  theory  as  being 
of  "zeroth"  order  in  dispersion  and  nonlinearity.  The  dispersive  index  M 
is  the  power  of  (h2^  2/3x2)  above  that  of  long  wave  theory.  The  nonlinear 
index  N  is  the  power  of  n/h0  above  that  of  long  wave  theory.  Thus  M=0 
corresponds  to  u1/2,  N=0  to  s1 ,  M=N=0  to  eui/2.  Tote  that  for  situations 
where  nonlinear  and  dispersive  effects  are  nearly  balanced,  this  theory 
extends  one  order  beyond  the  Boussinesq  Equations  (in  each  of  M,  !J).  For 
linear  dispersive  waves  this  theory  (when  -th  order  :oef f icients  from  2. I5 
are  used)  is  three  orders  beyond  the  3oussinesq  Equations.  This  theory 
recovers  long  wave  theory  when  M=0,  unlike  som^  versions  of  the  Boussinesq 


equations. 
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Fifnire  3  -  Dispersion  relations  Tor  linear  water  waves 


form  and  approximate  derivatives  with  respect  to  space  and  tine  by  centra! 
differences . 
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With  subscripts  denoting  the  spatial  position,  and  superscripts  the  tine 
level,  the  simplest  difference  from  of  (3.1)  that  places  both  n  and  q  at  the 
same  grid  points  is 

n+1  T7i—  1  n  n 

nn  -  nn  Vl  '  V’ 

- +  h  — d— i - =  r 

2At  o  2Ax 

1-3  o', 


TH+ - 


n-1 


2  At 


n+1 


-  n 


n-1 


2  At 


Upon  substitution  of  a  linear  wave  *  exp  i(uJt-kx)  for  n  and  q  we  derive: 

sin2(wAt)  ,  sin2(kAx)  ,  x 

- rAtT2-  =  gno  [A FT  (3-3 

Equation  (3.3)  is  neutrally  stable  for  j At | <j Ax!  (ghc)-1 /2,  in  that  w  is 
real  for  all  real  k.  The  dispersion  relation  for  At=Ax  (gh0)“1/z  is 

ui2/k 2  =  gh  (3.  M 

o 


which  is  exact  for  arbitrary  wave  numbers  and  numerical  resolution  Ax.  '■Then 
!  At | *|  Ax|  (ghQ) /2,  we  can  derive  the  relation 


U>2 

P5T 


gh  ( AO2;  +  0[  ’  Ax/h 
o'  o 


This  equation  contains  second  order  numerical  dispersion  unless 
|  At  I  =(gho  )~i /2[  Ax|  .  It  thus  appears  that  the  choice  of  '  At!  =(gho)_1  /  Ax; 
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will  likely  minimize  numerical  dispersion  for  fairly  long  waves.  ror  long 
linear  waves  this  choice  lies  at  the  boundary  of  instability,  however,  and 
we  must  examine  stability  in  a  little  detail  for  the  more  general  problem. 
To  do  this,  let  us  consider  the  effects  of  nonlinearity  and  dispersion 
separately . 

Consider  dispersion  first.  The  second  order  linear  set  of  equations 
takes  the  dimensionless  form: 


3j]  _3u  1_  33u 
3t  3x  15  3x^ 


2_  [,7  _  2  jfu,  +  3n 
3t  5  3x2 1  3x 


“i 


(  *5 


Distances  are  measured  in  units  of  h  and  soeeds  in  units  / er.  •  Again, 

o  o 

letting  the  field  he  represented  as  proportional  to  exp  (i(oJt-kx)  ]  ,  we  can 
derive  the  following: 

o 

1  +  p—  P 

sin2(wAt)  =  (~~)2  sin2  (kAx)  { - F— - }  (3.'T'' 

AX  1  +  $  P 

where  ?  =  [l-cos(kAx) j/( Ax)2  >  0.  Consequently,  the  expression  within  the 
brackets  of  (3.7)  never  exceeds  unity.  For  (At/ Ax' 2  <1,  the  PHS  of  (3.^) 
does  not  exceed  unity,  so  a  is  real,  i.e.  we  again  have  neutral  stability 
for  all  k ,  Ax. 

After  expanding  (3.7),  one  can  derive  the  numerical  dispersion  relation 
(f  or  !  At  |  =!  Ax |  ) : 


u)2  =  k2 


1  + 


A_  k2 

15 


1  +  4  k2 


+  k“*(Ax)2  +  higher  order  terms 


(Ax)2 
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The  first  term  of  the  right  han-i  side  of  (3.3.1  gives  the  exact  dispersion 
relation  for  (3.6).  The  lowest  order  numerical  dispersion  errors  are 
proportional  +  o  k"*(  Ax)  2,  with  a  small  coefficient,  in  contrast  to  a  choice 
|  At  <!  Ax|  ,  vr.ich  would  introduce  ar.  error  0[k2( Ax)2!  with  a  not  so  sraal  1 
coefficient  unless  At  is  almost  Ax. 

Other  formulations  of  the  problem.,  besides  being  of  lower  order, 
sometimes  choose  a  velocity  variable  in  such  a  way  that  'accidentally) 
constrains  practical  solutions.  ”cr  example,  if  qc  'the  same  us  uQ  in 
the  low  order  formulation)  is  taken  to  be  the  dependent  velocity  variable, 
the  linearised  Roussinesq  equations  are  still  in  the  form  of 

% 


(3.6),  but  with  u  replaced  by  q  and  different  numerical  coefficients:  ir. 


!3.6),  -1/13  goes  to  +1/3  and  2/5  gees  to  zero.  Then  (;.~'  becomes 


sin2  ( oiAt ) 


r  At  \  2 
'  Ax  ' 


sin2'/.  Ax' 


Here  the  factor  within  the  brackets  approaches  -=°  as  (Ax)2mE,  and  so 
sin2(oiAt)  <0  and  u  can  have  an  imaginary  part,  presumably  unstable.  Only  if 
At  goes  tc  cero  at  least  as  fast  as  'Ax)J  is  stability  assured.  This  makes 
numerical  work  very  expensive  for  equations  having  more  dispersion  in  forms 
like  the  first  of  (3.6)  than  in  forms  like  the  second,  because  more  time 
steps  must  be  taken.  The  Korteveg-deVries  Equation  has  expensive 
dispersion,  and  typical  finite  difference  solutions  even  as  coarse  as 
Ax=0.25  require  time  steps  the  order  of  IT-4  (See  Vliegenthart  l?'’!  for 
fuller  discussion).  The  Regularized  long  Wave  Equation,  on  the  other  hand, 
has  its  dispersion  in  a  form  akin  to  the  second  of  (3.6),  and  numerical  work 
can  be  done  at  At=Ax  (e.g.  3or.a,  et  al  19?0). 
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The  inclusion  of  some  nonlinear  terms  to  (3.o)  for  simple  waves  changes 

the  stability  analysis  only  by  a  little  when  the  waves  are  long.  Then  one 

can  approximate  the  nonlinear  terms  as  the  product  of  a  slowly  varying  **i®li 

and  a  potentially  unstable  rapidly  varying  field.  The  slowly  varying  field 

is  incorporated  in  the  long  wavelength  wave  speed,  resulting  in 

o 


sin2(uiAt)  =  (Cy~' )  sin2  1:  As  {- 
Ax 


LS  ‘ 


+  -  P 

0 


} 


where  c  =  ,’gh}1/2  ±  u  includes  nonlinear  contributions  and  can  exceed  unity. 
Thus  (cAt/Ax)2  is  a  factor  that  can  exceed  unity  when  At  =  Ax.  'lote,  however 


1+  -Ir  P 


that  this  does  not  imply  instability,  for  sin2  (kAx)  { - ^ — }  is  bounded 


1+  —  P 


0 

by  a  number  less  than  unity,  the  bound  being  dependent  on  (Ax)2.  For 
example,  for  Ax=l/°,  typical  of  the  computations  run  to  date,  this  factor 
never  exceeds  0.2.  Hence,  c  can  exceed  unity  by  a  comfortable  margin 
without  necessarily  giving  instabilities.  The  stability  analysis  given  here 
is  not  complete,  however.  Computational  instability  has  arisen  for  some 
cases  involving  very  high  waves.  The  conditions  under  whi^h  instability 
occurs  will  be  identified  later. 

B.  Boundary  Conditions 

In  all  calculations  run  to  date,  rigid  impermeable  side  and  end  walls 
bound  the  fluid.  For  Ia.place's  Equation  (2.7'  to  possess  valid  solutions  we 
must  specify  a  local  property  of  b  everywhere  along  a  closed  boundary.  Here 
the  boundary  conditions  at  the  lower  boundary  are  specified  by  the  form  of 


(2.9)«  The  top  boundary  condition  is  specified  by  (2.5)  which  says 


in  essence,  that  the  pressure  is  constant  there.  At  the  ends  the  physics 
demands  that  the  normal  components  fluid  velocity  match  the  wall  velocity 
a*  each  point.  This  can  be  transformed  to  conditions  on  n  and  q  demanded 
by  this  formulation.  Hiking  walls  to  be  vertical  and  moving  horizontal ly , 


u  (wall,v,t) 


U  ( t ' 


1 2 


where  U(t)  is  the  wall  speed. 

We  further  need  a  condition  on  n  at  the  wall.  This  is  most  easily 
derived  in  a  Lagrangiar.  formulation.  One  equation  of  motion  is: 

32x(a,b,t)  3x  ,  32y  (a,h,t )  i  3y  1  3p(a,b,t )  _  n  , , 
3t 2  3a  '  St2  ^  3a  p  3a 


where  a  and  b  are  the  initial  coordinates  of  the  fluid  particle  now  at 
(x,y).  Let  us  evaluate  the  terms  for  the  fluid  surface  particle  at  the 
wall.  When  the  wall  is  nonaccelerating  the  first  term  vanishes.  When  the 
fluid  surface  is  horizontal  at  the  initial  time,  the  third  term  vanishes. 
Apart  from  the  possibility  that  the  fluid  at  the  surface  accelerates  at 
exactly  -g,  this  demands  that  3y/3a=C,  which  is  the  same  as: 

(3.13) 

3x 

We  know  that  ug  (wall,  t)  =  U,  Using  (3.11) 

q  =  U  (3.1!*) 

5 

at  the  wall,  so  long  as  it  is  nonaccelerating.  We  do  not  consider 
accelerating  end  walls  in  this  report. 
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\  Initial  Conditions  and  Time  Levels 

All  calculations  are  started  by  specifying  the  wall  geometry  and 
values  of  n(x,t  =  O',  q  (x,t  =  •?)  throughout  the  region.  From  these  one 

5 

calculates  u  from  f 2.13^  and  u  from  (2.14)  at  t=C.  Values  of  n  and  o  at 

•s 

t=i(l/2)A*  are  formed  by  one-sided  time  differencing  and  the  fluxes  entering 
(2. If'  and  (2.1l'  at  t=+(l/2)  At  are  found.  The  fields  at  t=-(i/2'At  and 
fluxes  at  t=+(l/2)  At  are  used  to  advance  the  calculation  to  t=(3/2)At. 
"luxes  at  t=(3/2)At  and  fields  at  t=(l/2)At  are  used  to  advance  to  t=(3/2'At 
ar.d  the  process  is  repeated  indefinitely.  For  some  runs  the  process  is 
continued  to  t  =  (>;+(l/2)  )At,  ar.d  then  the  time  interval  is  reversed,  i.e.  , 
dt+-dt,  using  the  fields  at  t  =  (  ?.'+l/2)  At  as  initial  conditions.  Ihe  tine 
levels  or.  the  return  toward  t=0  are  t  — 1At ,  ar.d  lie  mi  dvay  between  the 
advancing-*, ime  levels  until  t=?  is  reached. 

There  appears  to  be  some  advantage  in  the  procedure  of  initially  moving 
to  two  adjacent  time  levels  symmetrically .  Other  procedures  such  as  a 
forward  or  backward  whole  time  step  difference  gave  some  alternating  grid 
point  noise  which,  of  course,  the  pure  leapfrog  method  does  not  suppress  as 
time  advances.  Trapezoidal  corrections  to  the  leapfrog  method  suppress  this 
oscillation,  but  also  damp  waves  a  little  and  decrease  the  accuracy  of 
solitary  waves.  I  have  seen  no  evidence  of  statistically  significant 
alternating  grid  point  noise  ir.  the  the  fields  when  the  calculations  are 
started  symmetrically,  as  described  above. 
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4 


where  a.  is  an  initial  amplitude,  and  a  is  related  to  it  by: 


"Hhe  relation  between  n  and  q  in  ^4.1)  and  between  a  and  a-  in  (U»2)  holds 
for  ‘he  (lowest  order'1  analysis  of  a  solitary  wave  found  in  Lar.b  (1932). 

The  computational  box  contained  13U0  grid  points;  time  was  advanced  through 
Ilf? I*  time  steps.  Tee  step  sices  were  Ax  =  At  =  l/°,  so  that  the  hex  length 
is  approximately  1=230  and  time  advanced  to  approximately  llQ.  feme  r.r,- 
solitary  wave  disturbance  was  present  1  r.  the  initial  conditions  for  all 
runs,  hut  the  solitary  wave  soon  outran  other  disturbances. 

Various  diagnostics  were  performed  to  determine  the  properties  ;f  this 
wave.  The  location  and  value  of  the  crest  elevation  was  found  by  fitting  a 
parabola  through  the  highest  elevation  and  its  neighbors,  from  which  a  wave 
speed  was  determined.  Various  Integral  properties  of  the  solitary  wave, 

such  as  its  mass,  of”  ndx,  vere  determined  by  numerical  integration 

n 

crest 

from  the  crest  to  the  right  side  of  the  box.  The  two  linear  conserved 
quantities  in  (2. A)  and  (2.5),  mass  and  velocity,  are  found  to  be  conserved 
properly  in  the  calculations. 

Specifically,  from  (2.1)  the  conserved  quantity  is  the  total  area  above 
still  water  level;  this  cemes  from  multiplying  (2.1^  by  b  and  integrating 
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Because  J=0  at  the  end  vails,  x  .  .  .  an.-}  x,  ^  .  only  the  first  tern  of 

right  left 

(4.3)  is  nonzero.  Integration  over  tine  gives 

x^.  h 

J  '  (nb)dx  =  Const 
'  left 


The  integral  of  (4.4)  is  monitored  at  each  step  of  the  calculation,  and  is 
found  to  be  constant  (to  the  four  significant  figures  printed  out  in  every 
calculation).  From  (2.5)  the  conserved  quantity  is  an  integral  of  q  , 


i.e. 


3  fright  ,  l/3 H.  2 

q  dx  +  |gn  -  -(-ry)2 
'left  3  2  3.  x 
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where  in  (4.5)  we  have  dropped  terms  involving  q_  and  n'  at  the  ends. 


because  they  vanish  there.  Ihe  integral  over  -“ime  produces: 
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Equation  (4.6)  says  that  at  each  tine  the  total  "velocity”  in  the  box 
/q  ix  varies  only  by  fluxes  entering  from  the  end  walls  (the  right  hand 
side  of  (4.6)).  Before  distrubances  reach  the  walls  we  find  that  fa  dx 
is  conserved  to  four  significant  figures.  .After  n*0  at  an  end  wall,  we  find 
the  differences  between  the  two  sides  of  (4.6)  are  snail  and  can  be 
explained  as  roundoff  and  truncation  errors. 

Figure  4  shows  the  wave  speed,  figure  5  the  total  solitary  wave  nass, 
and  figure  6  the  total  energy  of  the  solitary  wave.  Ihe  square.'  represent 
runs  with  second  order  dispersion  (see  2.16);  the  triangles  represent  third 
order  and  fourth  order  dispersion  runs  (2. I7  and  C.i?).  Because  these  are 
solitary  waves,  where  nonlinear  effects  balance  dispersive  ones,  the  formal 
accuracy  of  the  analysis  on  which  the  computations  are  based  is  limited  to 
second  order,  however,  and  improving  the  dispersive  contribution  al~ne 
should  not  effect  the  formal  accuracy  of  the  solutions.  From  "igs 
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is  clear  that  the  numerical  calculations  are  accurate  -  more  so  than  "secor. 


order"  solitary  wave  theories.  The  second  order  dispersion  gives  the  most 
accurate  results,  likely  a  fortuitous  result,  and  one  which  I  do  not 
emphasize.  'Jo  matter  how  the  dispersion  is  treated,  the  solitary  waves  are 
accurate.  The  highest  run  using  second  order  dispersion  shown  in  the 
figures  occurs  when  in  (U.l)  a. =0.5,  giving  a  wave  that  asymptotically  has 
an  amplitude  of  0.5^53.  When  a-=C.6  was  attempted,  the  calculations 
became  unstable.  Using  fourth  order  dispersion,  a; =0.6  leads  to  a  wave 
having  an  amplitude  of  0.72k.  The  run  for  a=0.7  became  unstable.  ”o 
attempts  at  finding  out  where  runs  with  third  order  dispersion  became 
unstable  were  made. 

The  resolution  Ax=At=l/8  is  sufficient  to  make  truncation  errors  small 
Results  at  resolution  Ax=At=l/k  show  slight  departures  from  those  displayed 
Runs  at  Ax=At=l/l6  are  unstable  at  lower  amplitudes  than  at  Ax=At=l,/8,.  Thi 
result  may  be  caused  by  roundoff  error,  and  might  not  hold  for  calculations 
carried  out  with  more  precision  than  the  6-7  places  used  'single  precision) 
Where  results  of  runs  using  Ax=At=l/R  and  Ax=At=l/l6  can  be  compared,  they 
agree  within  roundoff  errors.  This  is  in  accordance  with  the  linear 
analysis  leading  to  (3.2),  for  (1/36)  k4  (Ax!'2  is  very  small,  even  if 
the  estimate  of  appropriate  values  to  use  for  k4  are  not  so  small.  If  we 
assume  that  truncation  errors  are  proportional  to  (Ax'2,  the  small 
differences  observed  between  a  run  with  Ax=l/1  and  Ax=l/°  indicate  that 
truncation  differences  between  runs  with  Ax=l/S  ar.d  Ax=l/l6  should  be 
masked  by  roundoff  error,  as  observed. 
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expansion:;  are  from  lon^uet  -Higgins  a  mi  Fenton,  l'>74.  "Tie  first  order  wave 
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expansions  are  longuet, -lliggi  ns  and  Fenton,  197*'.  "Vie  first  order  energy 


Txamnles 


5. 

While  it  is  not  difficult  to  run  numerous  examples  of  water  waves  that 
have  different  boundary  geometries  and  initial  conditions,  this  paper  limits 
discussion  to  four  different  types  that  span  a  variety  of  conditions.  In 
subsection  A  an  example  of  solitary  waves  colliding  head-on  is  shown. 

In  B  examples  of  essentially  dispersive  wave  trains  are  shown.  In  C 
examples  of  solitary  waves  propgating  in  a  converging/diverging  channel  are 
shown.  In  D  an  example  of  long  waves  of  depression  coupled  with  an  unduiar 
bore  is  shown.  These  examples  are  designed  to  illustrate  the  capability  cf 
the  numerical  calculations  to  treat  waves  for  which  dispersive  effects  are 
larger  than,  less  than,  and  comparable  to  nonlinear  effects. 

A.  Solitary  Wave  Collision 

The  past  twenty  years  has  witnessed  a  surge  of  interest  in  whether  the 
solitary  waves  described  by  various  model  equations  are  solitons.  cere  I 
define  solitons  as  solitary  waves  that  emerge  from  a  collision  with  each 
another  having  the  identical  structure  that  they  had  prior  to  colliding  fan 
individual  soliton  may  be  retarded  or  advanced  during  the  collision).  The 
two  first  order  theories  that  have  been  used  to  study  colliding  solitary 
waves  in  water  are  the  Korteweg-deVries  equation 


3n  |  3n 
3t  3x 


=  o 


and  the  regularized  long  wave  equation 


3n  an  3_  _3n  1  33n 

3t  3x  +  2  3x  3x  z3t 


(5.2' 


The  formal  accuracies  of  (5.1'  and  (5.2)  are  the  same.  Ihe  first 
are  C (ep*/2),  t,he  next  is  and  the  last  is  C(eu^/^).  At 


two  terms 
2 (  eu1 /2' 
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the  solution  is  3/3t  =  -  3/3x,  so  that  in  the  last  term  of'  (  =  .  1 )  -3/3t  may 
be  substituted  for  3,/3x  to  give  (5.2)  without  modifying  the  formal  accuracy. 
Yet  the  solitary  waves  of  (5.1)  are  solitons  (Zahusky  and  Hruskal 
while  those  of  (5.2^  are  not  quite  solitons  'Hona,  et  a  1.  19s',  fart a r~ 1 1  i 
19"8,  Lewis  and  Tjon  1979). 

Here  I  describe  calculations  designed  to  see  whether  solitary  waves  in 
the  higher  order  theory  are  solitons.  The  still  water  depth  and  the  breadth 
of  the  channel  are  constant.  The  initial  "orm  of  each  wave  is  chosen  *0 
have  a  waveform  that  is  solitary -wave-like,  i.e. 

n  =  a  sech2a(x-x  )  (-.'' 

max 


where  x  marks  the  initial  location  of  a  wave  crest,  and  a  and  a  are 
adjustable  parameters,  as  in  Eq.  (4.1).  1ft er  some  experimentation,  I  nave 
found  that  the  choice  of  a  and  a  that  produces  acceptably  small  disturbances 
other  than  a  propagating  solitary  wave  differs  a  little  from  that  of  (L,?'. 
Here  a  solitary  wave  is  chosen  having  one  of  the  amplitudes  listed  in 
Table  5  Longs et -Higgins  and  Fenton  (197L),  'or  which  selected,  essentially 

exact,  solitary  wave  properties  are  known.  Among  these  are  the  speed  7, 

00  00 

f  ^  r  9 

the  mass  jn  dx,  and  the  potential  energy  ^  j  dx.  These  latter  two  are 

“00  — oo 

sufficient  to  determine  a  and  a  for  the  particular  solitary  wave  in  mind. 

The  initial  value  of  qg  is  found  by  setting  3/3t  =  ±F  3/3x  in  '2.c-''  'the  - 
sign  for  an  intended  right-going  wave  and  the  +  sign  for  left -goir. v ’  ,  and 
solving  the  resulting  quadratic  equation  for  qg  =  q3'F,n,n'). 
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Figure  'T  shows  plots  of  the  elevation  above  still  water  of  a 
disturbance  which  started  with  a  pair  of  solitary -like  waves.  The  wave 
initially  at  x  =  50  approaches  an  amplitude  of  0.1t"0  prior  to  the 
collision.  Hie  wave  initially  at  x  =  150  approaches  an  amplitude  of  O.38ol 
before  the  collision.  During  the  collision  the  elevation  reaches  D.635D, 
somewhat  more  than  the  sum  of  the  amplitudes  of  the  colliding  waves.  3y  the 
end  of  the  experiment  the  rightward  wave  had  reached  an  amplitude  of  0.1?<52 
(and  was  still  rising  slowly),  and  the  leftward  wave  had  reached  an 
amplitude  of  0.3851  (and  was  also  still  rising). 

Figure  8  takes  the  same  data  as  used  to  plot  ^ig.  "r,  expands  vertical 
scales  by  a  factor  of  3  and  clips  off  the  wave  crests.  Although  some 
oscillatory  wave  trains  are  visible  that  have  their  origin  at  the  initial 
condition,  the  major  oscillatory  wave  trains  that  fill  the  region 

-x. 

between  the  solitary  waves  at  t  >50  emanate  from  the  collision.  These 
waves  have  some  energy  that  the  solitary  waves  are  leaving  behind. 

Consequently,  the  solitary  waves  cannot  be  solitons  (strictly  speaking,  no 

more  than  one  of  them  can  be). 

The  question  remains  whether  the  existence  of  the  oscillatory  wave 
trains  that  show  the  solitary  waves  to  be  not  exact ly  solitons  can  be  an 
artifact  of  the  finite  difference  numerical  treatment.  One  potential  source 
of  error  could  he  what  numerical  analysists  refer  to  as  truncation  error, 
which  tends  to  be  proportional  to  (Ax)2  for  small  Ax.  To  rule  out  this 
possibility,  the  calculations  were  run  with  Ax  =  At  =  1/4,  i.e.  at  one  half 
the  resolution  of  Figs.  7-8.  The  oscillations  between  the  solitary  waves 
had  about  the  same  amplitudes  and  phases  at  both  numerical  resolutions. 
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For  example,  the  first  few  minima  and  maxima  to  the  right  of  the  larger 
solitary  wave  at  the  top  of  ^igs.  have  the  values  -T.CThp,  +9.0C9<:, 

-0.0C63,  +0.273+,  -O.Ohuv,  +0.001"’.  The  coarser  resolution  run  gives  the 
corresponding  values  of  -0.C0S'\  +  0.012'7,  -0.007®,  +0.2033,  -C.00<-'?,  0.0CT1. 
While  it  is  true  that  the  coarser  resolution  oscillations  are  a  little  more 
intense  than  the  finer  resolution  ones,  they  are  not  the  factor  1  expected 
under  the  hypothesis  that  the  existence  of  the  oscillatory  train  is  due  to 
truncation  error  (the  reason  the  oscillations  are  a  little  stronger  is 
likely  due  to  the  fact  that  with  the  sane  initial  conditions  the  solitary 
waves  are  a  little  higher  before  collision  for  Ax  =  At  =  1/1,  than  for 


Ax  =  At  =  1  /  ? ) 

To  see  whether  the  treatment  of  dispersive  terms  would  influence  the 
propagation,  a  run  at  Ax  =  At  =  1/3  with  fourth  order  dispersion  (not  the 
second  order  for  Figs.  7-8)  was  made.  Again,  the  oscillatory  wave  train  is 
present,  with  almost  exactly  the  same  amplitudes  and  nearly  the  same 
location  of  maxima  and  minima  as  in  Figs.  'T~8. 

Finally,  a  test  was  made  to  see  whether  any  diffusive  artifacts  were 
introduced  into  the  computer  program.  Figure  9  shows  the  results  of  a 
calculation  run  to  a  time  of  100,  then  time  reversed.  This  involves  a  total 
of  loOO  time  steps.  (As  for  Fig.  8,  the  amplitudes  are  clipped).  To  the 
scale  of  the  figure  the  calculations  are  perfectly  reversible.  The 
uppermost  profile  shown  i.n  Fig.  9  would  he  identical  to  the  initial  profile 
if  the  calculations  were  perfectly  reversible.  In  fact,  these  two  differ  by 
j(10-14)  throughout  the  computational  region.  We  ascribe  these  small 
differences  to  roundoff  error  in  running  through  loCO  time  steps  (SCO 
forward  and  800  back)  with  single-precision  accuracy  ( 6-~  significant 
figures).  The  phases  of  the  crests  ar;  remarkably  will  preserved.  For 
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example,  the  crest  of  the  leftward  propagating  wave,  which  started  at  x  = 
150.000  returned  to  x=  150.015,  after  the  round  trip  travel  distance  of  23^. 
We  may  conclude  that  these  colliding  solitary  waves  are  almost  solitons,  hut 
not  quite. 

3.  Dispersive  Wave  Trains 

The  wave  model  derived  here  is  designed  to  produce  solutions  for  linear 
and  not  so  linear  'wave  trains  that  are  not  necessarily  lor.ff.  To  see  how 
faithfully  the  model  represents  such  wave  trains,  consider  the  disturbance 
produced  from  the  following  initial  condition: 

n  =  0.08  (x-25)  sech2  0.08{x-25) 
q  =  n 

The  nondimensionalization  is  again  the  obvious  h  =  g  =  1.  This  tvte 
of  disturbance  resembles  that  of  an  impulse  distributed  over  a  spatial 
region  of  the  order  of  unity.  Figure  10  shows  the  results  of  the 
calculations  running  for  a  non-iimensional  time  of  100  and  time  reversed  for 
the  following  time  interval  of  100  (1600  time  steps  in  all).  The  result  is 
a  dispersive  wave  train.  The  longer  waves  travel  faster  than  the  shorter 
ones,  and  so  these  longer  waves  appear  at  the  front  of  the  train.  As  in  all 
of  the  figures  in  this  paper,  the  time  and  space  scales  are  set  so  that  a 
feature  of  a  disturbance  {a  crest,  for  example^  lines  up  with  a  slope  of 
exactly  unity  if  it  travels  with  a  speed  of  unity  (in  dimensionless  units; 
unity  corresponds  to  the  long  wavelength  limit  of  linear  wave  propagation, 

/gh^  in  iimensional  units).  As  Fig.  10  shows,  the  longer  waves  at  the 

front  of  the  train  have  speeds  slightly  less  than  unity,  while  the  shorter 
waves  travel  slower.  Waves  near  the  rear  of  the  train  are  very  short, 
having  length  scales  cons i lerably  less  than  the  depth.  (The  best  measure  of 
a  length  scale  of  a  periodic  wave  is  its  inverse  wave  number  ,  rather 
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Figure  10  -  Waves  from  an  impulse-1 ihe  initial  waveform,  with  time  reversal. 
The  time  scale  is  the  same  as  the  space  scale,  so  that  the  time  interval 
between  profiles  is  6.25  and  speeds  of  ±1  would  show  up  along  lines  oriented 
at  ±U5°.  Time  advances  to  t=TOO,  and  receeis  to  t=2C0.  The  elevation 
scale  is  0.05  per  interval:  the  initial  waveform  has  a  maximum  of  a  little 
less  than  0.05.  The  calculations  are  run  with  second  order  dispersion  ar.d 
Ax=At=l/8. 
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than  its  wavelength  A.  This  makes  the  measure  of  the  smallness  of  the  wave 
scale  equal  to  kh  ,  a  number  that  is  unity  when  the  waves  are  intermediate 
between  being  shallow  water  waves  and  deep  water  waves.  In  Fig.  10  ‘he 
waves  near  the  hack  of  the  forward  going  train  are  essentially  deep  water 
waves.)  Thus,  the  wave  model  sees  essentially  deep  water  waves,  even  though 
its  design  involves  an  expansion  in  a  parameter  that  vanishes  only  in  the 
limit  of  shallow  water  waves.  How  well  the  model  treats  deep  water  waves  is 
still  a  matter  of  conjecture,  because  sufficient  testing  has  not  yet  beer, 
performed.  It  is  comforting,  however,  that  we  see  the  deep  water  waves,  and 
the  dispersion  relations  shown  in  Fig.  3  say  that  waves  should  be  pretty 
well  represented  out  to  kh  as  large  as  0  or  more  (these  are  really  deep 
water  waves). 

As  in  the  calculations  that  produce  Fig.  9,  the  code  is  reversible, 
apart  from  round-off  error.  The  topmost  profile  is  identical  t  the  initial 
lowermost  profile  to  within  a  few  parts  in  IT1’.  The  maximum  elevation  of 
the  topmost  profile  occurs  at  the  location  x=  20. 9b,<,  where  it  was  initially 
at  25.965,  thus  indicating  very  good  phase  stability. 

In  order  to  make  more  quantitative  statements  about  how  well  the  model 
handles  dispersive  wave  trains,  calculations  v»re  performed  with  the 
following  initial  profile: 

n  =  0.5  sech2  0.2(x-12.5)  cos  2(x-12.*'  ('5.5' 

and  the  initial  surface  velocity  set  from  using  (2.5^,  with  F  being  the 
anticipated  wave  speed  of  a  periodic  wave  having  wave  number  2.0.  The 
resulting  disturbance  is  almost  entirely  right-going,  as  expected.  Figure 
11  shows  33  profiles,  equally  spaced  in  time,  running  from  t  =  0  to  t  =  c0. 
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~igare  11  -  Dispersive  wave  packet.  The  tine  scale  is  the  same  as  t 
scale,  so  that  the  time  interval  between  profiles  is  l.;t25.  The  el 
scale  is  0.15  per  interval;  the  intital  waveform  has  a  maximum  of  0. 
Tabulations  are  run  with  fourth  order  dispersion  and  Ax=£t=l/°. 
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Both  an  individual  wave  crest  and  -he  entire  disturbance  travel  at  speeds 
less  than  unity.  The  crests  travel  faster  tha-  the  group,  entering  at  the 
left  and  disappearing  at  the  right.  This  is  a  graphic  demonstration  that 
for  water  waves  the  phases  travel  faster  than  the  group. 

The  computational  data  that  make  up  Fig.  11  were  examined  in  some 
detail.  To  round-off  error,  the  largest  wave  crest  travels  at  the  linear 
speed  of  waves  of  wave  number  2.0.  Again  to  round-off  error,  the  entire 
group  travels  at  the  group  speed  of  linear  waves.  This  latter  result  is  no 
surprising,  even  though  the  wave  train  is  so  short  that  the  notion  of  a 
group  speed  (in  the  linear  theory)  is  a  little  fuzzy.  The  former  result, 
that  the  phase  of  the  largest  wave  traveled  at  the  linear  phase  speed, 
suprises  me  a  little.  Because  the  waves  near  the  center  of  the  train  are 
not  small-amplitude,  one  might  expect  that  they  would  move  with  speeds  a 
little  larger  than  the  predictions  of  linear  theory.  "or  example,  the  cres 
of  a  nonlinear  periodic  wave  described  by  Stokes  wave  theory  should  travel 
approximately  2  percent  faster  than  linear  waves  of  the  same  wavelength  for 
the  example  shown.  The  precision  of  the  calculations  of  crest  speed  is 
about  0.5  percent.  Why  the  largest  crest  travels  closer  to  what  linear 
theory  predicts,  rather  than  what  nonlinear  theory  predicts,  is  an 
unanswered  question  at  this  time.  It  may  have  to  do  with  the  shortness  of 
the  train,  so  that  the  largest  crest  is  intimately  connected  with  small 
amplitude  waves;  it  may  have  to  do  with  defining  a  zero  velocity,  which 
differs  in  the  initial/boundary  problems  solved  here  from  the  steady-state 
problem  solved  by  Stokes  and  others. 
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As  a  test  of  one  aspect  of  the  model,  the  data  used  in  generating 
Fig.  11  were  examined  a  little  more  carefully  to  see  whether  one  feature  of 
high,  not-so-long  waves  that  is  obvious  to  any  observer  of  waves  —  their 
sharp-crestedness  —  becomes  more  apparent  as  the  wave  amplitude  increases. 
The  waves  of  Fig.  11  are  more  nearly  deep-water  waves  than  shallow-water 
waves,  and  their  slope  is  the  better  measure  of  nonlinearity  than  their 
amplitude-depth  ratio.  Here  the  maximum  slope  is  about  0.3,  and  so  the 
waves  are  not  "small  amplitude";  they  are  also  not  "near-breaking",  which 
would  call  for  maximum  wave  slopes  of  the  order  of  0.6  (tan  30°).  Figure 
12  is  a  blowup  of  the  data  from  the  uppermost  profile  of  Fig.  II  (at  the 
-op)  and  of  the  central  profile  (at  t  =  25).  When  large  (the  center  of  the 
figure),  the  waves  are  somewhat  sharp-crested.  From  this  ve  can  conclude 
chat  the  wave  model  reproduces  the  sharp-crested  feature  possessed  by  high 
water  waves,  even  when  the  wavelength  is  not  so  long.  Whether  all  features 
cf  high,  not-so-long  waves  are  reproduced  from  this  long-wave  model  is  not 
yet  possible  to  say. 

C.  Waves  in  channels  of  variable  breadth 

The  wave  model  incorporates  a  variable  channel  breadth,  under  the 
circumstance  that  a  measure  of  the  length  over  which  the  breadth  changes 
substantially  is  nuch  greater  than  the  horizontal  scale  of  the  wave.  A 
computational  channel  was  set  up  to  replicate  conditions  of  experiments  by 
Char.g,  et  al.  (1979) •  The  channel  is  10  times  wider  or  narrower  at  one  end 
from  the  other.  We  connect  the  converging/diverging  parts  of  the  channel 
with  parallel  wall  geometry,  so  that  a  wave  can  travel  a  long  way  before 
encountering  the  variable  part  of  the  channel.  For  the  diverging  case, 
parallel  side  walls  occupy  the  region  between  x  =  9  and  x  =  12^.  The 
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profile  i.;  at  <.=50;  the  lower  profile 


walls  then,  diverge  from  :<  =  120  to  the  end  of  the  channel,  which  is  at  x  = 
200  (see  the  top  of  Fig.  1^1.  The  experiments  had  a  shorter  entrance 
sect! :n,  and  did  not  terminate  in  a  rigid  wall,  Before  the  wave  bounces  off 
the  end  wall,  however,  the  computations  should  mirror  the  experiments. 

They  do.  Figure  13  displays  the  propagation  of  a  wave  that  is 
designed  to  be  a  solitary  wave  initially,  as  in  earlier  displays.  The 
initial  conditions  are  set  to  be  the  same  as  those  of  the  leftward- 
propagating  wave  of  Fig.  10  (except  that  it  is  right -going ' ,  i.e.,  with  an 
amplitude  that  stayed  constant  near  0.336  at  the  wave  approached  the 
°r.  trance  to  the  diverging  part  of  the  channel  at  x  =  120.  Tine  wave 
generates  a  significant  reflected  wave  as  it  first  passes  through  the 
1i verging  section  of  the  channel.  It  is  here  that  reflections  are  expected 
to  he  the  largest,  for  the  nondimens  Iona l  length  associated  with  changes  of 
channel  breadth  are  the  largest  near  the  entrance  section.  The  solitary 
wave  gets  smaller,  and  acquires  an  oscillatory  tail.  It  also  slews  down; 
this  is  illustrated  more  clearly  in  Fig.  Ik,  which  uses  the  same  data  as  for 
Fig.  13,  hut  views  the  scene  from  the  coordinate  system  which  moves  with 

the  speed  of  the  long  wavelength  limit  of  linear  wave,  /qh  .  The  wave 
travels  considerably  faster  than  this  limiting  speed  at  first,  and  then 

slows  down  as  it  becomes  weaker,  still  exceeding  /gh ^ .  The  propagation 

speed  of  the  wave  crest  lies  close  to  that  of  a  solitary  wave  of  the  same 
amplitude.  Figure  15  expands  the  view  over  the  section  in  which  Chang,  et. 
al. 's  experiments  were  performed.  They  used  wave  probes  to  give  a  time 
history  of  elevation  at  one  of  several  -elected  locations.  They  also 
compare  their  observations  with  a  ’'crtewem-ieVries  ‘■heory.  f * i r.v  >r.  i 
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Figure  lU  -  Solitary  wave  propagating  to  and  through  a 
The  data  are  the  sane  as  for  fig.  13.  hut  are  viewed  f 


system  moving  to  the  right  with  a  dimensionless  speed 
coordinate  systems  being  aligned  at  t=0. 


McDonald  ( 1 9 8 2 )  compare  their  experiments  with  theories.  Let  it  suffice  to 
say  that  the  theory  given  here  reproduces  the  essential  features  of  the 
experiment,  such  as  the  os  illatory  tail  behind  the  wave  crest  generated  as 
the  wave  passes  through  the  diverging  section,  while  the  Korteweg-deVries 
theory  does  not. 

Figure  16  replicates  conditions  of  another  experiment  of  Chang,  et  al. 
(1979),  this  example  showing  waves  in  a  converging  channel.  Because  the 
physical  still-vater-depth  happened  to  he  different  in  this  'converging 
channel)  experiment  from  the  previous  (diverging  channel)  one,  the 
dimensionless  geometry  is  different.  Here  the  channel  entrance  section 
extends  to  x  =  lHO,  and  the  10-fold  converging  part  occupies  only  the  region 
between  x  =  lUO  and  x  =  200.  The  solitary  wave  ge^.s  bigger  'and  travels 
faster)as  it  propagates  in  the  converging  section.  Dramatic  reflected  waves 
and  oscillatory  wave  trains  following  the  major  dost rubar.ee  a:e  absent  from, 
this  scene.  There  is  a  long  positive  tail  to  the  wave,  however,  which  is 
consistent  with  the  experiments,  but  inconsistent  with  Korteweg-deVries 
theory  (see  Witting  ar.d  McDonald,  19?2). 

Two  asymptotic  theories  speak  to  the  behavior  of  waves  in  channels  of 
gradually  varying  breadth:  the  first,  derived  by  Creen  { 13  3  ?  \  predicts 
that  the  amplitude  of  the  disturbance  a  is  related  to  the  channel  breadth  h 
by: 

a  -  V1/2  '  = 

The  second,  discussed  in  detail  by  Miles  (l?’T9',  who  independently  derives 
earlier  results,  gives: 

a  a  -a-2/  3  ’  -  .  ' 

The  first  result  assumes  that  the  waves  can  be  described  by  linear  long  wave 
t.oeory,  and  that  the  disturbance  is  very  narrow  compare  J.  to  the  distance 
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DIMENSIONLESS  DISTANCE 


Figure  16  -  Solitary  wave  propagating  through  a  converging  channel.  The 
wave  was  launched  at  t  =  0,  x  =  25,  with  an  amplitude  of  0.2.  The  elevation 
scale  is  0.15  per  interval.  Calculations  are  rur.  with  second  dispersion  and 
Ax  =  At  =1/8. 


51 


1 


over  which  Che  channel  varies  substantially.  The  second  result  assumes  that 
the  waves  can  be  described  by  a  Ko rt ewe g— deVries  theory,  and  that  the 
channel  changes  breadth  sc  gradually  that  not  only  is  there  room  for  the 
wave  itself,  but  also  for  the  ’wave  to  continually  reform  itself  so  as  to 
remain  a  solitary  wave.  The  wave  propagations  shown  in  Figures  13-16,  and 
also  from  other  geometries  and  initial  conditions,  give  results  that 
generally  lie  between  the  extremes  of  Fq.  5.5  -  5*6,  and  cast  doubts  as  to 
whether  in  practice  the  requirements  of  the  asymptotic  theories  apply. 
Briefly,  it  appears  that  extremely  long  channels  ' L  =  hundreds)  may  be 
required  to  allow  sufficiently  gradual  changes  in  breadth.  7or  cases  run  in 
moderate-length-channels  (L  =  200)  there  is  not  enough  room.  Moreover,  in 
situations  where  a  channel  has  a  converging/diverging  section  with  entrance 
and  exit  sections  identical,  we  find  irreversible  behavior,  i.e.  ,  the 
emerging  waves  are  smaller  when  they  leave  the  variable  section  chan  they 
were  when  they  entered.  The  calculations  show  that  this  is  principally  the 
result  the  of  reflected  waves  that  escape  to  infinity  whenever  a  wave  passes 
through  a  section  of  variable  breadth.  3oth  asymptotic  theories  incorrectly 
predict  reversibility,  at  least  in  the  sense  that  the  amplitude  of  the 
emerging  wave  should  be  identical  to  that  of  the  entering  wave. 

D.  Undular  bores  and  long  waves 

Computational  experiments  were  run  to  see  whether  it  would  be  possible 
to  simulate  the  laboratory  experiments  of  Favre  (1935),  which,  even  today, 
form  the  most  complete  set.  One  aspect  was  not  modelled  —  the  manner  by 
which  the  bores  were  generted.  Favre  pushed  water  into  one  end  of  his  box; 
here  a  "dam"  is  broken,  the  water  filling  the  computational  region  being 


initially  still, 


Figure  17  simulates  one  of  Favre* s  experiments 


DIMENSIONLESS  DISTANCE 

Figure  IT  -  Tndular  bore  and  long  wave  disturbance  from  the  breaking  of  a 
gentle-faced  dam.  The  time  scale  is  the  same  as  the  space  scale,  sc  that 
the  time  interval  between  profiles  is  11. 3o25.  The  elevation  scale  is  0.; 
per  interval;  the  initial  elevation  at  the  left  is  0.1.  The  calculations 
are  run  with  fourth  order  dispersion  and  Ax  =  At  =  0.3125. 


(see  his  Fig.  Ul).  The  water  at  the  left  has  a  dimensionless  height  of  0.1 
Two  waves  are  launched  from  the  initial  condition:  a  long  wave  that  travel: 
to  the  left,  unsteepens,  and  reflects  from  the  wall  at  x  =  0.  This  wave 
could  Just  as  well  he  described  by  long  wave  theory  as  by  anything  more 
sophisticated.  The  initially  long  wave  travelling  to  the  right,  however, 
steepens  for  a  while,  and  then  starts  to  generate  a  wavetrain.  The  entire 
structure  to  the  right  of  x  =  ill  is  a  bore  that  soon  becomes  unduiar.  but 
never  reaches  a  steady  state.  The  maximum  elevation  of  the  undulations 
becomes  about  0.1  at  late  times. 

The  experiments  of  Favre  show  the  undulations  to  form  at  much  earlier 
times,  so  calculations  were  performed  with  an  initially  sharper  gradient  of 
elevation.  Figure  18  shows  the  results.  The  unsteepening  wave  of 
depression  travelling  to  the  left  is  now  accompanied  by  an  oscillatory  wave 
train,  generated,  presumably,  by  the  initial  sharp  gradient.  The  unduiar 
bore  is  qualitatively  the  same,  but  the  undulations  begin  to  form  earlier. 
Indeed,  the  number  of  waves  in  the  train  at  any  time  closely  remsembles  the 
number  in  the  wave  probe  records  shown  by  Favre  at  locations  comparable  in 
space-time. 

One  interesting  feature  of  the  calculations,  also  seen  by  Witting 
(1975),  is  that  the  water  depth  at  the  location  of  the  initial  disturbance 
rapidly  goes  to  the  value  given  by  long  wave  theory,  until  reflections  from 
end  walls  intrude.  For  the  conditions  of  Figure  17  this  is  not  surprising, 
for  the  wave  is  at  least  initially  long.  For  the  conditions  of  Figure  1^, 
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rur.  with  fourth  order  dispersion  and  Ax  =  At 


assumptions  of  the  the  long  wave  theory  are  violated  even  initially,  and  so 
the  nar.ner  by  whin',  the  system  adjusts  to  the  long  wave  results  is  an 
interesting  problem.  I  have  no  idea  how  to  go  about  deriving  this  observe  1 
result  from  analysis. 

To  see  the  undular  bore  structure  a  little  better,  the  data  of  Fig.  IS 
are  replotted  from  the  frame  of  reference  moving  to  the  right  with  a  speed 

of  /ghQ.  Plotted  are  the  initial  profile,  and  every  fourth  profile  of 

Fig.  IS.  This  corresponds  to  times  of  U6,  93,  139,  135,  291  and  2'78.  The 
qualitative  features  correspond  closely  to  calculations  using  the  Korteveg- 
deVries  Equation  (see  Peregrine,  1366  and  Vliegenthart ,  197l).  '’or  example, 
the  line  connecting  particular  crests  near  the  front  is  nearly  straight,  its 
slope  decreasing  as  time  increases.  Favre’s  experiments  at  early  times  show 
the  same  behavior.  At  late  times,  however,  Favre’s  data  show  "saturation", 
the  lead  few  waves  having  fche  same  amplitude,  only  the  waves  further  back 
shoving  the  evolving  sloping  line.  Thus,  the  results  of  the  computations 
giver,  here  do  not  model  this  feature  of  the  experiment  any  better  than  do 
other  theories. 

In  an  attempt  to  bette-  model  the  experiments,  some  dissipation  was 
added.  The  form  of  model  dissipation  is  that  of  the  trapezoidal  correction 
of  the  basic  leapfrog  method,  performed  at  each  time  step.  Figure  TO  shows 
the  results  for  the  sane  conditions  as  in  T!g.  1?.  ^igure  20  looks  like  a 
damped  picture  of  Fig.  19,  i.e.  ,  the  waves  are  smaller,  but  the  overall 
envelopes  of  the  profiles  are  qualitatively  the  same.  It  is  true  that  there 
are  differences  in  the  model  dissipation  and  the  laboratory  dissipation. 
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Figure  1?  -  Detailed  view  of  an  undular  here.  The  data  are  samples  of  the 
same  data  shown  in  Fig.  18,  but  are  viewed  from  the  coordinate  system  moving 
to  the  right  with  a  dimensionless  speed  of  unity,  the  coordinate  systems 
being  aligned  at  t=0.  Shown  are  the  profiles  at  every  118th  time  step,  and, 
from  bottom  to  top,  occur  near  t=0,  16,  03,  139,  18,5,  031,  and  D'’?. 
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figure  20  -  Detailed  view  of  an  undular  bore  with  dissipation.  Initial  data 


and  running  conditions  are  the  same  as  for  Figs.  17-1?.  The  only 


differences  between  this  and  Fig.  1?  are  due  to  the  effects  of  running  with 
trapezoidal  correction  at  each  time  step. 
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Nevertheless,  it  is  likely  that  the  difference  between  experiment  and  theory 
is  due  to  something  other  than  dissipation,  perhaps  surface  tension.  This 
contrasts  with  the  case  treated  in  the  last  subsection,  where  the 
experiments  definitely  favor  this  theory  over  a  Korteweg-deVries  theory. 

6.  Discussion 

The  end  result  of  the  work  reported  here  is  a  new  model  of  water  waves 
that  can  describe  a  wide  variety  of  propagtion  situations  accurately  and 
efficiently.  The  following  factors  are  responsible  for  this  development: 

1.  The  model  uses  exact  prognostic  equations  in  conservation  form, 

Bqs.  (2.1)  and  (2.5).  \  very  general  derivation  of  (2.5)  is  given  in 

Witting  and  McDonald  (19'32).  To  ny  knowledge  it  has  not  been  used  in  wave 
modeling  before  now. 

2.  Higher  order  expansions  than  used  before  connect  the  velocity 
variables  that  appear  in  the  governing  equations.  This  allows  the  model  a) 
to  incorporate  long  wave  theory  exactly,  b)  to  give  a  good  representation  of 
waves  all  the  way  out  to  khQ  exceeding  8,  thus  including  both  shallow 
water  and  essentially  deep  water  waves  in  the  sane  model,  and  c)  to 
represent  fairly -long  nonlinear  waves  to  one  order  better  than  Boussinesq. 

3.  The  model  employs  a  numerical  method,  i.e,  pure  leapfrog,  that 
gives  no  unwanted  numerical  diffusion.  The  time-stepping  procedures  are 
simple  enough  to  analyze  in  some  detail  and  to  implement  efficiently  on 
vector  computers. 
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U.  The  model  Takes  a  time  step  equal  to  a  space  step  (in 
nondimensional  units  for  which  the  linear  long  wave  speed  is  unity).  This 
allows  for  efficient  machine  computations,  unlike  methods  developed  for  The 
Korteweg-deVries  Equation.  Moreover,  this  procedure  removes  any  spurious 
numerical  dispersion  at  order  k2(Ax)2. 

5.  Finally,  the  diagnostic  equations  are  cast  in  a  form  such  that  only 
tri^iiagonal  matrix  equations  need  to  be  solved.  A  very  fast,  fully 
vectorized  algorithm  is  then  used  to  invert  the  matrices. 

Running  times  for  the  waves  model  on  the  TI^ASC-7  computer  are 
approxintely  20  msec  per  time  step  for  a  computational  region  containing 
1600  grid  points.  A  run  to  t  =  150  with  resolution  Ax  =  At  =  l/°  takes 
about  2h  sec.  About  half  of  the  running  time  is  spent  in  collecting 
diagnostic  information  and  can  be  eliminated,  if  desired.  Ir.  summary , 
computer  costs  to  nun  the  model  are  very  small. 

?!o  detailed  comparisons  have  been  made  with  other  numerical  models  for 
specific  problems.  In  general,  though,  it  is  clear  from  Figs.  (2-6)  that 
the  model  is  considerably  more  accurate  than  Boussinesq  models,  and  can  deal 
with  higher  and  shorter  waves.  The  flexibility  of  the  model  to  treat 
different  physical  conditions  is  illustrated  by  Figs.  (7-20)  with 
accompanying  text.  Consequently,  the  model  represents  a  great  improvement 
over  other  expansion-type  models.  Fxcept  for  breaking  waves  ard  other 
extreme  situations  that  it  cannot  reach,  the  unified  waves  model  can  treat 
the  same  problems  that  non-expansion-tvpe  models,  such  as  lagrangian  and 
Marker-in-cell  models,  can  treat,  at  a  tiny  fraction  of  the  cost.  In 
practice,  these  more  complicated  models  cannot  employ  a  resolution  fine 
enough  to  compete  with  the  new  waves  model. 
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